Formation energy and interaction of point defects in two-dimensional colloidal crystals 
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The manipulation of individual colloidal particles using optical tweezers has allowed vacancies 
to be created in two-dimensional (2d) colloidal crystals, with unprecedented possibility of real-time 
monitoring the dynamics of such defects (Nature 413, 147 (2001)). In this Letter, we employ molec- 
ular dynamics (MD) simulations to calculate the formation energy of single defects and the binding 
energy between pairs of defects in a 2d colloidal crystal. In the light of our results, experimental 
observations of vacancies could be explained and then compared to simulation results for the in- 
terstitial defects. We see a remarkable similarity between our results for a 2d colloidal crystal and 
the 2d Wigner crystal (Phys. Rev. Lett. 86, 492 (2001)). The results show that the formation 
energy to create a single interstitial is 12% — 28% lower than that of the vacancy. Because the pair 
binding energies of the defects are strongly attractive for short distances, the ground state should 
correspond to bound pairs with the interstitial bound pairs being the most probable. 



I. INTRODUCTION 



A revival of general interest in point defects in two- 
dimensional (2d) systems has been fueled by demonstra- 
tion that vacancies can be created in 2d colloidal crystals 
through manipulation of individual particles with optical 
tweezers^. Because of the large size of the colloidal par- 
ticles, the structural and dynamical properties of these 
point defects can be monitored with video microscopy 3 . 
Interest in point defects in solids is widespread, includ- 
ing defects in ordinary materials as well as in quantum 
crystals such as He and Wigner crystals. In quantum 
crystals the point defects are believed to occur at finite 
concentrations at any nonzero temperature. There is also 
speculation that even at zero temperature point defects 
can exist, while at higher concentrations they may lead to 
a supersolid phase^£. The role of point defects in melt- 
ing in a two-dimensional system has been investigated 
theoretically^!*^, but an experimental investigation is 
usually hampered by the low concentrations of defects, 
unless the temperature is close to the melting points. 

In this work, we report on the first accurate Molecu- 
lar Dynamics (MD) calculations of the formation energy 
of single point defects and interaction of pair point de- 
fects in a 2d screened Coulomb interaction colloidal sys- 
tem. The system is formed by colloidal particles, taken 
as identical spheres of radius a, suspended in a solvent 
(generally water) and confined between two parallel solid 
surfaces. When immersed in this solvent, the colloids ac- 
quire a large charge Z due to dissociation of endgroups 
from their surface. Counterions thus generated ensure 
charge neutrality, forming a cloud around each colloid 
that makes the Coulomb interaction shorter. The col- 
loids are free to move in 2d and interact through pair- 
wise Yukawa type potentiat^-i 3 ". The Hamiltonian for 
this system is 
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where the first term in the right is the kinetic energies, 
the second is the screened Coulomb colloid-colloid in- 
teraction and the third term is the interaction between 
the colloids and the neutralizing background of positive 
charges 14 Ub = — 27r&A/ag, where b = (2/v / 3) 1 ^ 2 ao is the 
lattice space and ao is the average separation between 
colloids (this is defined for the triangular lattice with 
vector translation (6,0) and (6/2, by/3/2) with a = 1/^/p 
where p is the 2d colloid number density), A is the screen- 
ing length, e is the dielectric constant of the medium 
and Z\ = Z*f(a/X), where Z* is the normalized charge 
and f(x) = sinh(x)/a; is a function that describes the 
effect of the nonzero radius a of the colloidal particles. 
The energy, length, temperature and time are in units of 
Eq = (Z\e) 2 /ea, where a = 1.1/im (typical lattice space 
for experimental systems), To = Eq/kb {kb Boltzmann 
constant) and to = (Eo/ma 2 ) -1 ^ 2 , respectively. 



II. COMPUTATIONAL DETAILS 

The system with point defects (vacancy or interstitial) 
is modeled by removing (adding) a colloidal particle from 
(to) the most stable 2d lattice, which is the triangular lat- 
tice. We place the defect at the center of the simulation 
box to avoid complications arising from lattice relaxation: 
a single 6 coordinated vacancy at the central lattice site 
or a threefold centered interstitial colloid in one of the 
central triangular unit cells. However, no constraint ex- 
ists which would restrict the center of each colloid to lie 
within its own Wigner-Seitz cell, i.e. to make the lattice 
relaxation locally. This means that the defects are free 
to move around and can change their symmetry during 
the system evolution in thermodynamic equilibrium. 



2 



In order to calculate the energy needed to create one 
single defect we perform two independent simulations at 
the same density and temperature: a simulation for the 
perfect and a simulation for the defective system. For 
the latter, after introducing the defect, we rescale the di- 
mensions of the simulation box by a factor / = \Jnd/n p 
to reset the system to the original density, with rid and 
n p being the number of colloids for the defective and the 
perfect system, respectively. This is performed to cir- 
cumvent the need of correcting the energy due to density 
change caused by inclusion of the defect. The difference 
between the energies of the defective and the perfect sys- 
tems is the energy needed to create the defect. Formally, 
we can define the number of defects Ndef as the number 
of colloids minus the number of lattice sites. Therefore, 
the formation energy of Ndef defects in a crystal with TV 
lattice sites is 
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FIG. 1: Formation energy of a vacancy and interstitial defect 
as a function of temperature at p = 0.954 (/im)" 2 . Computed 
for systems of 120 lattice sites. 



AE def = [E(N + Ndef) - E(N)](N + N def ), (2) 



where E(n) is the energy per colloid for a system con- 
taining n colloids. For a monovacancy or an interstitial 
defect, Ndef = N v = -1 and N de f = N mt = +1, respec- 
tively. 

MD calculations of the formation energy of point de- 
fects were performed for several system sizes n = 29, 30, 
31, 129, 130, 131, 269, 270, 271, 479, 480 and 481 colloids, 
which allows one to study finite-size effects at different 
densities. Larger system sizes would be unpractical com- 
putationally because high accuracy is needed in order to 
obtain the energy differences. 

We used colloids with radius a ~ 0.18 ^m, screening 
length k^ 1 ~ 0.39 /im, charge ~ 1650e and density vary- 
ing from p — 0.402 to 1.804 (/itm) -2 , corresponding to 
typical experimental data^. To give an idea of the ener- 
gies and temperatures involved in the calculations, for in- 
stance at p — 0.954(^m)~ 2 the units are E = 9.4 x 10~ 18 
Joules and T = 6.81 x 10 5 K. 

The initial positions for the colloids are the sites of 
a triangular lattice accommodated in a rectangular box 
with periodic boundary conditions to eliminate surface 
effects. The simulations were performed within the 
canonical ensemble keeping a constant system tempera- 
ture using the Berendsen's thermostat with coupling pa- 
rameter to an external bath tt ranging from 0.01 to 0.1—. 
The evolution of Newton's equation of motion is obtained 
with the four-order predictor-corrector algorithm. The 
time step varies from 2.5 x 10 -2 to 5.0 x 10~ 3 to since it 
has some scale dependence on the colloid density Ther- 
modynamic equilibrium was assumed to be achieved dur- 
ing the first 50000 time steps, after which the physical 
quantities were obtained by averaging over 700 blocks of 
10000 time steps. 



III. RESULTS AND DISCUSSION 

In the simulation, the formation energy of the vacancy 
and interstitial defects depends on the number of parti- 
cles, i.e. size of the simulation box. To eliminate the size 
dependence of the formation energy we extrapolate the 
simulation results to the thermodynamic limit by using 
the following formula 

AE^ f (p) = AEXef(p) + ^, (3) 

where the density-dependent parameters AE^^(p) and 
c(p) are determined by a linear least-squares fit to the 
MD calculations at different TV. The MD formation en- 
ergies for the point defects and x 2 of the fitted data as 
a function of p, at a fixed temperature and screening 
length, are shown in table I. The \ 2 of the fitted data is 
particularly good, indicating that the size dependence is 
well described by Eq. (3). The formation energies for in- 
terstitial defects are consistently lower, from 12% — 28%, 
than those for vacancies in the range of densities studied. 
Therefore, an interstitial defect is more stable and more 
likely to be created. 

The formation energy of the defects is practically 
temperature-independent in the temperature range 3. x 
10~ 5 T < T < 5. x 10~ 2 T below the melting point, as 
shown in Fig. 1, where the formation energy of the in- 
terstitials is smaller than that of the vacancies. However, 
when the kinetic energy starts to compete with the poten- 
tial energy close to the melting point (T > 5.0 x 10 _2 To), 
the formation energy of the defects becomes too noisy due 
to large fluctuations in the total energy and approaches 
zero within the statistical error. This suggests that sin- 
gle point defects may be easily created thermally in this 
temperature range and could play an important role in 
the melting mechanisms of 2d colloidal crystals. 
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TABLE I: Formation energy in reduced units obtained from MD simulations for a vacancy (left) and interstitial defect (right) 
for various system sizes denoted by N (number of lattice sites) and at various values of p in (/im)~ 2 . Quantities in parentheses 
are the estimated error in the last decimal place. Also given are the energies for an infinite system (AEy ac , AEf^ t ), c(p) and 
X 2 fitting parameters. The system temperature and the screening length are T = 3.74 x 10 -4 Tb and kT 1 — 0.39 fi , respectively. 



Vacancy 

p = 0.402 p = 0.589 p = 0.954 p = 1.804 

iV = 30 0.0261(5) 0.0543(4) 

N = 120 0.029(1) 0.0550(4) 0.1029(4) 0.1933(4) 

N = 270 0.029(2) 0.055(1) 0.1032(4) 0.1958(7) 

N = 480 0.029(3) 0.055(1) 0.1032(9) 0.1958(6) 

AE^ ac 0.029(1) 0.0551(4) 0.1033(4) 0.1969(7) 

c(p) -0.10(3) -0.02(2) -0.05(8) -0.4(1) 
X 2 0.086 0.032 0.023 0.6 



We now investigate the possible defect topologies as 
the system evolves in thermodynamic equilibrium, which 
is particularly relevant for defining the dynamics of the 
defects, as observed experimentally^ and recently corrob- 
orated by results from Brownian simulation methods^. 
According to our simulations, the initial configurations 
for the threefold centered interstitial and the sixfold va- 
cancy were found to relax into a configuration of lower 
symmetry, in agreement with simulation and experimen- 
tal observations. However, the aforementioned observa- 
tions concerning the dynamics of defects focus only on 
the symmetry and topology of the defect. In the follow- 
ing we discuss the dynamics of the defects in topological 
and energetic terms. 

As stated in Rcf. -, at finite temperatures this many 
body system vibrates around every local energy mini- 
mum due to thermal fluctuations. If the energy differ- 
ences between distinct local minima are small, the sys- 
tem can get enough energy to move to a nearby local 
minimum. As long as the system remains around a local 
energy minimum, the distortions in the lattice are elastic 
and the topological arrangement of the particles does not 
change. This allows us to calculate the system energy for 
each topology. For such a calculation, the defect must 
be tracked after the system has reached thermodynamic 
equilibrium. For that, we developed a code to perform 
a dynamical check of a list of neighbors of each colloid 
in a triangular lattice. This list is created (updated) 
by counting the sides of the polygons in the Voronoiii 
construction at each time step run. A defect is charac- 
terized by the presence of miscoordinated particles, i.e. 
particles whose number of neighbors is different from 6. 
Once the current topology of the defect is identified, the 
corresponding energy is recorded. We also calculated the 
time the defect remains in a given topology and the num- 
ber of transitions each defect performs between different 
topological configurations. 

Table II summarizes our findings for the dynamics of 
the defects. The topologies for the vacancy are: crushed 
vacancy (V2); symmetric vacancy (V3); split vacancy 
(SV); which were observed experimentally 1 , and an- 



Interstitial 

p = 0.402 p = 0.589 p = 0.954 p = 1.804 

N = 30 0.0238(5) 0.0430(4) 

N = 120 0.024(1) 0.0435(6) 0.0781(5) 0.1425(3) 

N = 270 0.024(2) 0.043(3) 0.0783(6) 0.1426(4) 

iV = 480 0.024(2) 0.043(6) 0.0785(5) 0.1428(5) 

AEZt 0.024(1) 0.0436(7) 0.0785(6) 0.1428(5) 

c(p) -0.00(3) -0.01(3) -0.0(1) -0.03(8) 

X 2 0.00056 0.046 0.018 0.049 



other one, a fourfold symmetric excited configuration 
(V 4 ), only observed recently^. For interstitial defects 
the topologies are: threefold symmetric interstitial (I3); 
twofold symmetric interstitial (I2); disjoint twofold sym- 
metric interstitial (hd); and a fourfold symmetric excited 
interstitial (I 4 ). The transition matrix (a stochastic ma- 
trix) in the upper part of Table II indicates that for a 
sufficiently long time, both vacancy and interstitial de- 
fects adopt the possible topologies many times. Each row 
of this matrix gives the probability of transition from 
a state (topology) to another one. As pointed out in 
RefJ^, the topologies V2, SV, I2 and l^d have a prefer- 
ential diffusion direction, with motion being a random 
walk along the main crystalline directions. In contrast, 
the V3, V 4 , I3 and J 4 have no preferential direction and 
may act in switching the direction of motion. Therefore, 
the diffusion process of our system, for a long time run, 
is isotropic. However, we observed in subsidiary simula- 
tions (results not shown here) that for short runs there is 
a very small number of transitions to topological config- 
urations responsible for changing the direction of motion 
(V3, V 4 , I3 and J 4 ), especially for the vacancy. As a re- 
sult, diffusion becomes one-dimensional, consistent with 
experimental data and simulation resultsi^. 

At the bottom of Table II, we show the times (t s ) for 
the defects in each topological configuration in equilib- 
rium, as well as the formation energy AE for each topo- 
logical configuration. The formation energy for the vari- 
ous topologies are very close to each other, being almost 
within the statistical error. However, one infers that the 
lowest formation energy for both vacancy and intersti- 
tial defects correspond to the configurations Sy and I3, 
where in average the defects remain most of the time (t s ). 
This result should be expected as these configurations 
are the most likely to be thermally activated since they 
need the smaller amount of energy to be created. On the 
other hand, the defect tends to spend a very short time 
in V4 and 7 4 configurations because they have the high- 
est formation energies. Note that l^d topology is just a 
transient variation of I2 topology. A direct comparison of 
the formation energies for each topological configuration 
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FIG. 2: Binding energy of a pair defect of vacancy-vacancy 
and interstitial-interstitial as a function of lattice spaces at 
p = 0.954 (pm)~ 2 . Computed for systems of 120 lattices 
sites. 

and its lifetime can lead to wrong interpretation of the 
results of Table II. Although the defects formation ener- 
gies between such topological configurations are similar, 
the average time of the defect in each topological config- 
uration is quite different due to differences in the energy 
barriers. For instance, though the formation energy of 
V2 is slightly smaller than that of V3 , the lifetime of V 2 is 
much larger. The reason is that the net transition proba- 
bility from V3 to SV is much larger than that from Vi to 
SV (about, 36%, see transition matrix), which indicates 
a smaller energy barrier between V3 and topologies 
making V3 less stable than Vz- One can make also the 
same discussion for the interstitials defect (right hand 
side of the table II) . 

Far from the melting point, an increase in tempera- 
ture causes the number of transitions to increase, includ- 
ing the low-frequency transitions between different topo- 
logical configurations. Close to the melting transition, 
the formation energy for different topological configura- 
tions is difficult to calculate owing to the large thermal 
fluctuations, with the differences in energy lying within 
statistical error. 

We also studied the interaction of a pair of defects as a 
function of the lattice separation. The difference between 
the energies to create two single defects separated by a 
given number of lattice constants and the energy to create 
two isolated single defects is the binding energy between 
two single defects, defined as follows 

E lnt = AE def {N def = 2) - 2AE def {N def = 1) (4) 

where the first term in the right AE de f(N de f = 2) 
corresponds to the energy to create two single defects 
separated by some lattice constants, while the last term 
2AE de f(N de f = 1) is the energy to create two isolated 
single defects. 



Fig. 2 shows the binding energy vs. lattice separation. 
In order to enforce accuracy in positioning the pairs of 
defects, we slow down the movement of the colloids by 
decreasing the temperature to ~ 10 -10 To. It was enough 
to hold the two defects in a fixed separation. Both the 
pairs of vacancy-vacancy and interstitial-interstitial are 
strongly attractive at short distances, with attraction go- 
ing to zero for distances greater than three lattice spaces. 
For short distances the binding energies for the defects 
are ~ —0.08 and ~ —0.04 for the vacancy and inter- 
stitial defect, respectively, being therefore higher than 
the value expected at the melting temperature at den- 
sity p — 0.954 (^m) -2 , which is ~ 0.01 according to Fig. 
1. We have also (not shown here) calculated the binding 
energy as a function of the density for just one lattice 
space separation between defects, which resulted always 
attractive in the range of densities p — 0.402 — 1.804 
(p,m)~ 2 . Since any attraction should suffice to permit 
recombination, our results suggest that the ground state 
energy of the 2d colloidal crystal may be dominated by 
pair binding of defects. Furthermore, as interstitial de- 
fects have the lowest excitation energy, we expect that 
the ground state of point defects should involve mainly 
interstitial pairs. Though the mechanisms of pair bind- 
ing of point defects in solids are still not fully understood, 
we note that our results are similar to those for intersti- 
tial defects in quantum crystals, such as vortex crystals 9 
and 2d Wigner crystals 8 , in which the defect pair is not 
sufficiently strong to yield a supersolid phase. For experi- 
mentalists working in the 2d colloidal system using video 
microscopy, we believe the results on the pair binding of 
defects provide an experimental challenge to observe for- 
mation of pair point defects near the 2d colloidal crystal 
melting. 

Before concluding, we comment on the entropic fac- 
tors, which are embedded in the calculations of the for- 
mation energies. Such factors could be neglected well 
below the melting temperature, as they are much less 
important than the energy terms. Close to the melt- 
ing temperature, entropic factors are likely to affect the 
numerical values of the formation energies, but the qual- 
itative features should be preserved. Indeed, this expec- 
tation appears to be fulfilled in the experimental system 
represented by the 2d colloidal crystal, since our simula- 
tions could explain the experimental observations. 

In conclusion, we have shown quantitatively the ef- 
fects from point defects (vacancy and interstitial) in 2d 
colloidal crystals, in which the energy to create a single 
interstitial defect is 12% — 28% lower than to create a sin- 
gle vacancy. The formation energies of these defects go to 
zero near the melting point, i.e. point defects can be eas- 
ily created thermally and should play a crucial role in the 
melting mechanisms for such crystals. We also confirmed 
previous results that the interstitial defects arc more mo- 
bile than vacancies, and provided an explanation based 
on energy calculations. Finally, we found that the inter- 
action between defects is strongly attractive, and conse- 
quently most defects will exist as bound pairs. We believe 
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TABLE II: The upper part is the normalized transition matrix for different topological configurations of the defects, during 
~ 7.0 x 10 MD steps after the system has reached equilibrium. In the bottom are displayed the time spent in each topological 
configuration (t s ) and the formation energy in reduced units of the defects, AE, for the 2d colloidal crystal at p = 0.954 
(/im)" 2 . Quantities in parentheses are the estimated errors in the last decimal place. The system temperature and screening 
length are T = 1.0 x 10~ 3 To and kT 1 = 0.39 /im, respectively. 







Vacancy 








Interstitial 






V 2 


V 3 


SV 


vi 




h 


1 2d 


h 


h 


v 2 


0.9987098 


0.0000154 


0.0012747 





h 


0.9977724 


0.0000890 


0.0021340 


0.0000046 


v 3 


0.0000175 


0.9995009 


0.0004666 


0.0000150 


hd 


0.0051196 


0.9947921 


0.0000883 





sv 


0.0009514 


0.0000255 


0.9990231 


0. 


h 


0.0007948 


0.0000007116 


0.9989618 


0.0002426 


vi 


0. 


0.0048446 





0.9951554 


h 


0.0000450 


0.0000064 


0.0054687 


0.9944798 


ts 


0.4068799 


0.0420110 


0.5509791 


0.0001300 


ts 


0.2618633 


0.0046109 


0.7024319 


0.0310939 


AE 0.1044(6) 


0.1046(6) 


0.1040(6) 


0.1086(6) 


AE 0.0798(6) 


0.0799(6) 


0.0790(6) 


0.0829(6) 



that our results may have important bearing on experi- 
mental works involving interfaces and solid surfaces. 
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